function f = custom_tpdf(x,v)
% custom_tpdf - PDF of Students t distribution
%
% FORMAT:       f = custom_tpdf(x, v)
%
% Input fields:
%
%       x           t-ordinates
%       v           degrees of freedom (v > 0, non-integer d.f. accepted)
%
% Output fields:
%
%       f           PDF of t distribution with v d.f. at points x
%
% Note: this function has been copied from the SPM2 package published
%       by the Wellcome Department, go here for details:
%       http://www.fil.ion.ucl.ac.uk/spm/

% Definition:
%-----------------------------------------------------------------------
% The Student's t-distribution with v degrees of freedom is defined for
% positive integer v and x in (-Inf,Inf), and has Probability Distribution
% Function (PDF) f(x) given by: (See Evans et al., Ch37)
%
%               gamma((v+1)/2)
%    f(x) = -----------------------  *  (1 + x^2/v) ^ -((v+1)/2
%           sqrt(pi*v) * gamma(v/2)
%
% This implementation is not restricted to whole (positive integer) df
% v, rather it will compute for any df v>0.
%
% Algorithm:
%-----------------------------------------------------------------------
% Direct computation using the beta function for
%       sqrt(pi)*gamma(v/2) / gamma((v+1)/2)  =  beta(v/2,1/2)
%
% References:
%-----------------------------------------------------------------------
% Evans M, Hastings N, Peacock B (1993)
%       "Statistical Distributions"
%        2nd Ed. Wiley, New York
%
% Abramowitz M, Stegun IA, (1964)
%       "Handbook of Mathematical Functions"
%        US Government Printing Office
%
% Press WH, Teukolsky SA, Vetterling AT, Flannery BP (1992)
%       "Numerical Recipes in C"
%        Cambridge

% Version:  v0.7f
% Build:    8110521
% Date:     Nov-05 2008, 9:00 PM CET
% Author:   Jochen Weber, SCAN Unit, Columbia University, NYC, NY, USA
% URL/Info: http://wiki.brainvoyager.com/BVQXtools

% argument check
if nargin < 2 || ...
   ~isnumeric(x) || ...
   ~isnumeric(v) || ...
    isempty(v) || ...
    any(isnan(x(:))) || ...
    any(isinf(v(:)) | isnan(v(:)) | v(:) < 0)
    error( ...
        'BVQXtools:BadArgument', ...
        'Missing or invalid argument given.' ...
    );
end
if isempty(x)
    f = [];
    return;
end
osize = size(x);
if numel(x) == 1 && ...
    numel(v) > 1
    osize = size(v);
end
if ~isa(x, 'double')
    x = double(x(:));
else
    x = x(:);
end
if ~isa(v, 'double')
    v = double(v(:));
else
    v = v(:);
end
if numel(v) ~= numel(x) && ...
    numel(v) ~= 1 && ...
    numel(x) ~= 1
    error( ...
        'BVQXtools:SizeMismatch', ...
        'Invalid number of elements in v.' ...
    );
end
if numel(x) ~= numel(v)
    if numel(x) == 1
        x = repmat(x, size(v));
    else
        v = repmat(v, size(x));
    end
end

% compute and reshape
f = ( (1 + x .^ 2 ./ v) .^ (- (v + 1) / 2) ) ./ ...
	 (sqrt(v) .* beta(v ./ 2, 1 / 2));
f = reshape(f, osize);
